DMRscaler defines differential methylation features iteratively,
expanding the size of the window for aggregating on at each step this
procedure allows a hierarchical structure to describe a region enriched
in differntially methylated CpGs. We look at our example region here
dmr_tree <- example_generate_dmr_tree(dmrscaler_result = dmrscaler_result,
layer=length(dmrscaler_result),
chr=example_region$chr,
start = example_region$start,
stop = example_region$stop)
diagonalNetwork(List = dmr_tree, fontSize = 12, fontFamily = "bold",
nodeStroke = "black", linkColour = "black", opacity = 1 )
## organize data used in tracks
group <- character(length=ncol(B))
group[cases] <- "case"
group[controls] <- "control"
delta_mean_B <- rowMeans(B[,cases]) - rowMeans(B[,controls])
which <- which(locs$chr==example_region$chr & locs$pos >= example_region$start & locs$pos <= example_region$stop)
Bsub <- B[which, ]
gr <- GRanges(seqnames = example_region$chr, ranges = IRanges(start = locs[which,"pos"], width = 1), mcols = Bsub )
## set up ideogram track
itrack<-IdeogramTrack(genome="hg19", chromosome = example_region$chr)
## set up genome axis track
gtrack<-GenomeAxisTrack()
## set up significance track
sig_track <- DataTrack(range=gr, data = -log10(locs$pval[which]), type=c("p"))
## set up beta value track
beta_track <- DataTrack(range = gr, data = t(Bsub), groups=group,name="Beta", type=c("a","g","confint"))
## set up diff beta track
diff_beta_track <- DataTrack(range = gr, data = delta_mean_B[which], type=c("a","g"))
## set up gene model track ###
genesub <- example_gene_anno_df(chr=example_region$chr, start=example_region$start, stop=example_region$stop)
if(!all(is.na(genesub))){
grtrack <- GeneRegionTrack(genesub, chromosome = example_region$chr, start = example_region$start, stop=example_region$stop, transcriptAnnotation="symbol" ,collapseTranscripts = "meta")
}else{grtrack<-GeneRegionTrack()}
## set up DMRscaler results layer tracks
layer_track <- list()
for(j in 1:length(dmrscaler_result)){
if(nrow(dmrscaler_result[[j]])==0){
layer_track[[j]] <- AnnotationTrack()
} else {
layer_track[[j]] <- AnnotationTrack(start = dmrscaler_result[[j]]$start,
end = dmrscaler_result[[j]]$stop,
chromosome = dmrscaler_result[[j]]$chr,
strand = "*", genome = "hg19", name=paste("layer", j, sep = "_"))
displayPars(layer_track[[j]]) <- list(cex.legend=0.7, cex.axis=0.7, cex.title=0.7, rotation.title=0, stackHeight = 1, shape="box")
}
}
plotTracks(c(list(itrack, gtrack, sig_track, beta_track, diff_beta_track, grtrack), layer_track), from = example_region$start-1, to=example_region$stop+1 )
